function [mpe,beta,delta] = estimate(model)

    % Set-up for state space
    glob.n          = 10;           % Number of nodes
    glob.nf         = 300;          % Number of points in fine grid
    glob.curv       = 0.3;          % Grid curvature for A on (0,1] (1 is no curvature)
    glob.spliorder  = 3;            % Order of splines
    glob.amin       = 0;            % Lower bound on A
    glob.amax       = 3e6;          % Upper bound on A

    % Model parameters: taken as given
    glob.T         = 60;           % Period in which agents die (80-20)
    glob.R         = 45;           % Period in which agents retire (65-20)
    glob.r         = 0.02;         % Risk-free rate
    
    if isequal(model, 'clno')
        glob.gammac    = 20000;        % Subsistence consumption
        glob.wt        = [ones(glob.R,1).*160656]./1633; % Avg posttax income is 160656 and hours are 1633
        glob.at        = 0.7*160656;   % Retirement income
        glob.tau       = 0;            % tax rate is zero since all is post-tax
        glob.min_age   = 20;           % Minimum age in the model
        glob.gammah    = 1880;         % Maximum hours
    elseif isequal(model, 'winner-pretax-1.5%')
        glob.gammac     = 20000;        % Subsistence consumption
        glob.wt         = [ones(glob.R,1).*194505]./1633; % Avg pre-tax income is 194505 (see Table 3) and hours are 1633
        glob.tau        = 0.3;         % tax rate, see figure A.6
        glob.at         = 0.7*194505;  % (Pre-tax) Retirement income
        glob.min_age   = 20;           % Minimum age in the model
        glob.gammah    = 1880;         % Maximum hours
    elseif isequal(model, 'winner-pretax-2.5%')
        glob.gammac     = 20000;        % Subsistence consumption
        glob.wt         = [ones(glob.R,1).*194505]./1633; % Avg pre-tax income is 194505 (see Table 3) and hours are 1633
        glob.tau        = 0.3;         % tax rate, see figure A.6
        glob.at         = 0.7*194505;  % (Pre-tax) Retirement income
        glob.min_age   = 20;           % Minimum age in the model
        glob.gammah    = 1880;         % Maximum hours
    elseif isequal(model, 'household-pretax-2.5%')
        glob.gammac     = 20000 * 2;    % Subsistence consumption (2-person HH)
        glob.wt         = [ones(glob.R,1).*381865]./(2 * 1633); % Avg pre-tax household income is 381865 (see Table 6) and hours (per person) are 1633
        glob.tau        = 0.3;          % implied tax rate, see Figure A.6
        glob.at         = 0.7*381865;   % (Pre-tax) Retirement income
        glob.min_age    = 20;           % Minimum age in the model
        glob.gammah     = 1880 * 2;     % Maximum hours (2-person HH)
    end

    % First guesses of model parameters: to be estimated
    param.beta     = 0.855;        % Consumption weight
    param.delta    = 0.015;        % Discount factor

    % Simulation parameters
    % For age distribution:
    age_data = csvread('./moments/age_distribution.csv',0,0);
    glob.age_bins       = age_data(:,1) - glob.min_age - 1;
    glob.age_probs      = age_data(:,2)./sum(age_data(:,2));
    glob.age_probs      = round(glob.age_probs,3);
    glob.Na             = 1000;
    glob.N              = round(sum(glob.age_probs)*glob.Na);         % Number of agents

    % For initial wealth:
    % Pre-lottery net wealth in the three age groups are 273081, 637677 and 900574
    glob.mean_wealth = sum(glob.age_probs(1:14))*273081 ...       % initial wealth
    + sum(glob.age_probs(15:34))*637677 + sum(glob.age_probs(35:end))*900574;

    % For prize distribution
    glob.prize_dist   = [.097 .815 .066 .015 .001 .006];          % weights on prizes
    glob.prizes       = [500 5000 60000 350000 ...              % possible prizes
    750000 1500000];

    % Weighting matrix and moments
    if isequal(model, 'clno')
        effects_data = [csvread('./moments/moments.csv',2,2); 40];
        glob.W       = diag(effects_data(1:end-1).^2);
        targets = csvread('./moments/moments.csv',2,1);
        glob.moments_data = [targets(:, 1)];
    elseif isequal(model, 'winner-pretax-1.5%')
        glob.W       = diag(ones(10,1));   % identity matrix
        targets = csvread('./moments/figure_1.csv',5,1);
        glob.moments_data = [targets(:, 1)];
    elseif isequal(model, 'winner-pretax-2.5%')
        glob.W       = diag(ones(10,1));   % identity matrix
        targets = csvread('./moments/figure_1.csv',5,1);
        glob.moments_data = [targets(:, 1)];
    elseif isequal(model, 'household-pretax-2.5%')
        glob.W       = diag(ones(10,1));   % identity matrix
        targets = csvread('./moments/figure_5.csv',3,1);
        glob.moments_data = [targets(:, 1)];
    end
    
    %% Set up grids/splines

    % State space for endogenous variable A
    Na              = glob.n;
    curv            = glob.curv;
    spliorder       = glob.spliorder;
    agrid           = nodeunif(Na,glob.amin.^curv,glob.amax.^curv).^(1/curv);  % Adds curvature

    % Function space and nodes (fspace adds knot points for cubic splines)
    fspace          = fundef({'spli',agrid,0,glob.spliorder});
    sgrid           = funnode(fspace);
    s               = gridmake(sgrid);
    Ns              = size(s,1);
    glob.fspace     = fspace;
    glob.s          = s;    
    glob.Ns         = Ns;

    %% Solve & simulate
    opts = optimset('fminsearch');
    opts.Display = 'iter';
    if isequal(model, 'clno')
        est_params      = zeros(2,1);
        est_params(1)   = param.beta;
        est_params(2)   = param.delta;
        m = @(e_p) min_distance(e_p,glob,param);
        opt_ests = fminsearch(m,est_params, opts);
    elseif isequal(model, 'winner-pretax-1.5%')
        est_params      = zeros(2,1);
        est_params(1)   = param.beta;
        est_params(2)   = 0.015;
        m = @(e_p) min_distance(e_p,glob,param);
        opt_ests = fminsearchbnd(m,est_params,[-Inf 0.015],[Inf 0.015], opts);
    elseif isequal(model, 'winner-pretax-2.5%')
        est_params      = zeros(2,1);
        est_params(1)   = param.beta;
        est_params(2)   = 0.025;
        m = @(e_p) min_distance(e_p,glob,param);
        opt_ests = fminsearchbnd(m,est_params,[-Inf 0.025],[Inf 0.025], opts);
    elseif isequal(model, 'household-pretax-2.5%')
        est_params      = zeros(2,1);
        est_params(1)   = param.beta;
        est_params(2)   = 0.025;
        m = @(e_p) min_distance(e_p,glob,param);
        opt_ests = fminsearchbnd(m,est_params,[-Inf 0.025],[Inf 0.025], opts);  
    end
    
    % estimated parameters
    beta = opt_ests(1);
    delta = opt_ests(2);
    
    %% Model fit
    
    % moments, paths, objective function at optimum
    [moments, ~] = solve_simulate(opt_ests,glob,param,-1,0);
    
    % figure
    plot(0:10, [0;moments])
    hold on;
    plot(0:10, [0;targets(:,1)])
    hold off;
    if isequal(model, 'clno')
        saveas(gcf, 'output/model-fit-clno.png')
    elseif isequal(model, 'winner-pretax-1.5%')
        saveas(gcf, 'output/model-fit-winner-pretax-1.5p.png')
    elseif isequal(model, 'winner-pretax-2.5%')
        saveas(gcf, 'output/model-fit-winner-pretax-2.5p.png') 
    elseif isequal(model, 'household-pretax-2.5%')
        saveas(gcf, 'output/model-fit-household-pretax-2.5p.png')   
    end
    
    %% Compute lifetime wealth effect at age 30, i.e. MPE
    
    [~, path]  = solve_simulate(opt_ests,glob,param,30 - 20,0);
    diff_earnings   = path.earn_win - path.earn_nowin;
    dL              = path.wealth_win(:,1) - path.wealth_nowin(:,1);
    dydL            = bsxfun(@rdivide, diff_earnings, dL);
    mpe             = mean(nansum(dydL, 2));
       
end

